Module 3: Introductory Analytics

library(tidyverse)
demographics = read.csv('data/demos_anonymized.csv')
ids = read.csv('data/ids_anonymized.csv')
model_variables = read.csv('data/model_variables_anonymized.csv')

Introduction

Regression

model_1 = lm(award_total_amount ~ libuser, 
             data = model_variables)

coef(model_1)
summary(model_1)

Note there are only two fitted values, one for library users and one for the non-users.

fitted(model_1)[1:10]  # first 10 predictions on current data set

Use modelr to add them to your data set.

library(modelr) # need to install?
model_variables %>% 
  add_predictions(model_1, var = 'Prediction') %>% 
  group_by(libuser) %>% 
  summarise(result = mean(Prediction))

The broom library also has functionality for many common models.

library(broom) # need to install?
augment(model_1)

Use ggeffects for a quick visual!

library(ggeffects) # need to install?
ggpredict(model_1) %>% 
  plot()

We can assess performance: - within a model (e.g. Adjusted \(R^2\)) - by comparing models (e.g. anova test) - with prediction on new data (e.g. RMSE)

model_2 = lm(award_total_amount ~ libuser + gender, 
             data = model_variables)

anova(model_1, model_2)

Other approaches (all must be done on same data!):

Adjusted R^2

adjr1 = summary(model_1)$adj
adjr2 = summary(model_2)$adj

round(c(adjr1, adjr2), 2)

AIC is uninterpretable on its own, but reflects the likelihood of the data given the model and complexity. When you have multiple AIC values, the lower value is the better model.

AIC(model_1, model_2)

Classification

model_logreg = glm(award_total_amount > 5e6 ~ age, 
                   data = model_variables,
                   family = binomial)

summary(model_logreg)

The caret package has a nice function to summarize classification performance via a variety of metrics.

It needs labeled factors to use though. Often you would have this anyway.

library(caret) # need to install?

predictions =  predict(model_logreg) > 0  # same as probability > .5
predictions = factor(predictions, labels = c('low', 'high'))

target = model_variables$award_total_amount > 5e6
target = factor(target, labels = c('low', 'high'))

confusionMatrix(predictions, target)

Python examples

Init

import pandas as pd
import numpy as np
import statsmodels


model_variables = pd.read_csv('data/model_variables_anonymized.csv')

Statsmodels

Statsmodels is Python’s attempt to do statistical modeling. It essentially follows R’s formula style.

import statsmodels.api as sm
import statsmodels.formula.api as smf

# specify the model, then fit
model_1 = smf.ols('award_total_amount ~ libuser', data=model_variables)
results_1 = model_1.fit()

# Inspect the results
print(results_1.summary())

model_2 = smf.ols('award_total_amount ~ libuser + gender', data=model_variables)
results_2 = model_2.fit()

# Inspect the results
print(results_2.summary())

print(sm.stats.anova_lm(results_1, results_2))

Classification

model_variables['award_high'] = 1*(model_variables['award_total_amount'] > 5e6) # 1* makes it a numeric

logit_mod = smf.logit('award_high ~ libuser + gender', data = model_variables)
logit_result = logit_mod.fit()

print(logit_result.summary())

logit_result.pred_table() / model_variables.shape[0]

LS0tCnRpdGxlOiAiTW9kdWxlIDM6IEludHJvZHVjdG9yeSBBbmFseXRpY3MiCm91dHB1dDogCiAgaHRtbF9ub3RlYm9vazogCiAgICBoaWdobGlnaHQ6IHB5Z21lbnRzCiAgICB0aGVtZTogc2FuZHN0b25lCiAgICBjc3M6IG90aGVyLmNzcwplZGl0b3Jfb3B0aW9uczogCiAgY2h1bmtfb3V0cHV0X3R5cGU6IGlubGluZQotLS0KCmBgYHtyIGluaXQsIGVjaG89RkFMU0V9CiMgdGhlc2Ugb3B0aW9ucyBhcmUgcHJpbWFyeSB1c2VmdWwgdG8gdGhlIGNyZWF0aW9uIG9mIHRoZSBodG1sIGRvY3VtZW50CmtuaXRyOjpvcHRzX2NodW5rJHNldCgKICBlY2hvPVQsIAogIGV2YWwgPSBGLAogIG1lc3NhZ2UgPSBGLCAKICB3YXJuaW5nID0gRiwgCiAgY29tbWVudCA9IE5BLAogIFIub3B0aW9ucz1saXN0KHdpZHRoPTEyMCksIAogIGNhY2hlLnJlYnVpbGQ9RiwgCiAgY2FjaGU9RiwKICBmaWcuYWxpZ249J2NlbnRlcicsIAogIGZpZy5hc3AgPSAuNywKICBkZXYgPSAnc3ZnJywgCiAgZGV2LmFyZ3M9bGlzdChiZyA9ICd0cmFuc3BhcmVudCcpCikKYGBgCgpgYGB7ciBjYXRjaHVwfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKZGVtb2dyYXBoaWNzID0gcmVhZC5jc3YoJ2RhdGEvZGVtb3NfYW5vbnltaXplZC5jc3YnKQppZHMgPSByZWFkLmNzdignZGF0YS9pZHNfYW5vbnltaXplZC5jc3YnKQptb2RlbF92YXJpYWJsZXMgPSByZWFkLmNzdignZGF0YS9tb2RlbF92YXJpYWJsZXNfYW5vbnltaXplZC5jc3YnKQpgYGAKCiMjIEludHJvZHVjdGlvbgoKCiMjIFJlZ3Jlc3Npb24KCmBgYHtyIG1vZGVsMX0KbW9kZWxfMSA9IGxtKGF3YXJkX3RvdGFsX2Ftb3VudCB+IGxpYnVzZXIsIAogICAgICAgICAgICAgZGF0YSA9IG1vZGVsX3ZhcmlhYmxlcykKCmNvZWYobW9kZWxfMSkKc3VtbWFyeShtb2RlbF8xKQpgYGAKCk5vdGUgdGhlcmUgYXJlIG9ubHkgdHdvIGZpdHRlZCB2YWx1ZXMsIG9uZSBmb3IgbGlicmFyeSB1c2VycyBhbmQgb25lIGZvciB0aGUgbm9uLXVzZXJzLgoKYGBge3IgZml0dGVkfQpmaXR0ZWQobW9kZWxfMSlbMToxMF0gICMgZmlyc3QgMTAgcHJlZGljdGlvbnMgb24gY3VycmVudCBkYXRhIHNldApgYGAKClVzZSBgbW9kZWxyYCB0byBhZGQgdGhlbSB0byB5b3VyIGRhdGEgc2V0LgoKYGBge3IgbW9kZWxyX2RlbW99CmxpYnJhcnkobW9kZWxyKSAjIG5lZWQgdG8gaW5zdGFsbD8KbW9kZWxfdmFyaWFibGVzICU+JSAKICBhZGRfcHJlZGljdGlvbnMobW9kZWxfMSwgdmFyID0gJ1ByZWRpY3Rpb24nKSAlPiUgCiAgZ3JvdXBfYnkobGlidXNlcikgJT4lIAogIHN1bW1hcmlzZShyZXN1bHQgPSBtZWFuKFByZWRpY3Rpb24pKQpgYGAKClRoZSBgYnJvb21gIGxpYnJhcnkgYWxzbyBoYXMgZnVuY3Rpb25hbGl0eSBmb3IgbWFueSBjb21tb24gbW9kZWxzLgoKYGBge3IgYnJvb21fZGVtb30KbGlicmFyeShicm9vbSkgIyBuZWVkIHRvIGluc3RhbGw/CmF1Z21lbnQobW9kZWxfMSkKYGBgCgpVc2UgZ2dlZmZlY3RzIGZvciBhIHF1aWNrIHZpc3VhbCEKCmBgYHtyIHBsb3RfcmVncmVzc2lvbn0KbGlicmFyeShnZ2VmZmVjdHMpICMgbmVlZCB0byBpbnN0YWxsPwpnZ3ByZWRpY3QobW9kZWxfMSkgJT4lIAogIHBsb3QoKQpgYGAKCgpXZSBjYW4gYXNzZXNzIHBlcmZvcm1hbmNlOgotIHdpdGhpbiBhIG1vZGVsIChlLmcuIEFkanVzdGVkICRSXjIkKQotIGJ5IGNvbXBhcmluZyBtb2RlbHMgKGUuZy4gYW5vdmEgdGVzdCkKLSB3aXRoIHByZWRpY3Rpb24gb24gbmV3IGRhdGEgKGUuZy4gUk1TRSkKCmBgYHtyIHJlZ19tb2RlbF9jb21wYXJpc29ufQptb2RlbF8yID0gbG0oYXdhcmRfdG90YWxfYW1vdW50IH4gbGlidXNlciArIGdlbmRlciwgCiAgICAgICAgICAgICBkYXRhID0gbW9kZWxfdmFyaWFibGVzKQoKYW5vdmEobW9kZWxfMSwgbW9kZWxfMikKYGBgCgpPdGhlciBhcHByb2FjaGVzIChhbGwgbXVzdCBiZSBkb25lIG9uIHNhbWUgZGF0YSEpOgoKLSBBZGp1c3RlZCBSXjIKCmBgYHtyIGFkanIyfQphZGpyMSA9IHN1bW1hcnkobW9kZWxfMSkkYWRqCmFkanIyID0gc3VtbWFyeShtb2RlbF8yKSRhZGoKCnJvdW5kKGMoYWRqcjEsIGFkanIyKSwgMikKYGBgCgoqKkFJQyoqIGlzIHVuaW50ZXJwcmV0YWJsZSBvbiBpdHMgb3duLCBidXQgcmVmbGVjdHMgdGhlICpsaWtlbGlob29kKiBvZiB0aGUgZGF0YSBnaXZlbiB0aGUgbW9kZWwgYW5kIGNvbXBsZXhpdHkuICBXaGVuIHlvdSBoYXZlIG11bHRpcGxlIEFJQyB2YWx1ZXMsIHRoZSBsb3dlciB2YWx1ZSBpcyB0aGUgYmV0dGVyIG1vZGVsLgoKYGBge3IgYWljfQpBSUMobW9kZWxfMSwgbW9kZWxfMikKYGBgCgoKIyMgQ2xhc3NpZmljYXRpb24KCmBgYHtyIG1vZGVsX2xvZ3JlZ30KbW9kZWxfbG9ncmVnID0gZ2xtKGF3YXJkX3RvdGFsX2Ftb3VudCA+IDVlNiB+IGFnZSwgCiAgICAgICAgICAgICAgICAgICBkYXRhID0gbW9kZWxfdmFyaWFibGVzLAogICAgICAgICAgICAgICAgICAgZmFtaWx5ID0gYmlub21pYWwpCgpzdW1tYXJ5KG1vZGVsX2xvZ3JlZykKYGBgCgoKVGhlIGNhcmV0IHBhY2thZ2UgaGFzIGEgbmljZSBmdW5jdGlvbiB0byBzdW1tYXJpemUgY2xhc3NpZmljYXRpb24gcGVyZm9ybWFuY2UgdmlhIGEgdmFyaWV0eSBvZiBtZXRyaWNzLgoKSXQgbmVlZHMgbGFiZWxlZCBmYWN0b3JzIHRvIHVzZSB0aG91Z2guICBPZnRlbiB5b3Ugd291bGQgaGF2ZSB0aGlzIGFueXdheS4KCmBgYHtyIGNhcmV0fQpsaWJyYXJ5KGNhcmV0KSAjIG5lZWQgdG8gaW5zdGFsbD8KCnByZWRpY3Rpb25zID0gIHByZWRpY3QobW9kZWxfbG9ncmVnKSA+IDAgICMgc2FtZSBhcyBwcm9iYWJpbGl0eSA+IC41CnByZWRpY3Rpb25zID0gZmFjdG9yKHByZWRpY3Rpb25zLCBsYWJlbHMgPSBjKCdsb3cnLCAnaGlnaCcpKQoKdGFyZ2V0ID0gbW9kZWxfdmFyaWFibGVzJGF3YXJkX3RvdGFsX2Ftb3VudCA+IDVlNgp0YXJnZXQgPSBmYWN0b3IodGFyZ2V0LCBsYWJlbHMgPSBjKCdsb3cnLCAnaGlnaCcpKQoKY29uZnVzaW9uTWF0cml4KHByZWRpY3Rpb25zLCB0YXJnZXQpCmBgYAoKIyMgUHl0aG9uIGV4YW1wbGVzCgojIyMgSW5pdAoKYGBge3B5dGhvbiBweV9pbml0LCBlbmdpbmUucGF0aD0gJy9Vc2Vycy9taWNsL2FuYWNvbmRhMy9iaW4vcHl0aG9uJ30KaW1wb3J0IHBhbmRhcyBhcyBwZAppbXBvcnQgbnVtcHkgYXMgbnAKaW1wb3J0IHN0YXRzbW9kZWxzCgoKbW9kZWxfdmFyaWFibGVzID0gcGQucmVhZF9jc3YoJ2RhdGEvbW9kZWxfdmFyaWFibGVzX2Fub255bWl6ZWQuY3N2JykKYGBgCgojIyMgU3RhdHNtb2RlbHMKCmBTdGF0c21vZGVsc2AgaXMgUHl0aG9uJ3MgYXR0ZW1wdCB0byBkbyBzdGF0aXN0aWNhbCBtb2RlbGluZy4gIEl0IGVzc2VudGlhbGx5IGZvbGxvd3MgUidzIGZvcm11bGEgc3R5bGUuCgoKYGBge3B5dGhvbiBzdGF0c21vZGVsc30KaW1wb3J0IHN0YXRzbW9kZWxzLmFwaSBhcyBzbQppbXBvcnQgc3RhdHNtb2RlbHMuZm9ybXVsYS5hcGkgYXMgc21mCgojIHNwZWNpZnkgdGhlIG1vZGVsLCB0aGVuIGZpdAptb2RlbF8xID0gc21mLm9scygnYXdhcmRfdG90YWxfYW1vdW50IH4gbGlidXNlcicsIGRhdGE9bW9kZWxfdmFyaWFibGVzKQpyZXN1bHRzXzEgPSBtb2RlbF8xLmZpdCgpCgojIEluc3BlY3QgdGhlIHJlc3VsdHMKcHJpbnQocmVzdWx0c18xLnN1bW1hcnkoKSkKYGBgCgpgYGB7cHl0aG9uIHJlZ21vZGVsXzJ9Cm1vZGVsXzIgPSBzbWYub2xzKCdhd2FyZF90b3RhbF9hbW91bnQgfiBsaWJ1c2VyICsgZ2VuZGVyJywgZGF0YT1tb2RlbF92YXJpYWJsZXMpCnJlc3VsdHNfMiA9IG1vZGVsXzIuZml0KCkKCiMgSW5zcGVjdCB0aGUgcmVzdWx0cwpwcmludChyZXN1bHRzXzIuc3VtbWFyeSgpKQpgYGAKCmBgYHtweXRob24gbW9kZWxfY29tcGFyZX0KcHJpbnQoc20uc3RhdHMuYW5vdmFfbG0ocmVzdWx0c18xLCByZXN1bHRzXzIpKQpgYGAKCgojIyMgQ2xhc3NpZmljYXRpb24KCgpgYGB7cHl0aG9uIGxvZ2lzdGljfQptb2RlbF92YXJpYWJsZXNbJ2F3YXJkX2hpZ2gnXSA9IDEqKG1vZGVsX3ZhcmlhYmxlc1snYXdhcmRfdG90YWxfYW1vdW50J10gPiA1ZTYpICMgMSogbWFrZXMgaXQgYSBudW1lcmljCgpsb2dpdF9tb2QgPSBzbWYubG9naXQoJ2F3YXJkX2hpZ2ggfiBsaWJ1c2VyICsgZ2VuZGVyJywgZGF0YSA9IG1vZGVsX3ZhcmlhYmxlcykKbG9naXRfcmVzdWx0ID0gbG9naXRfbW9kLmZpdCgpCgpwcmludChsb2dpdF9yZXN1bHQuc3VtbWFyeSgpKQpgYGAKCmBgYHtweXRob24gY29uZnVzaW9uX21hdHJpeH0KbG9naXRfcmVzdWx0LnByZWRfdGFibGUoKSAvIG1vZGVsX3ZhcmlhYmxlcy5zaGFwZVswXQpgYGAKCg==